1 - Número de SADs preditas não refutadas

Figura 1.1 logito(número de SADs não refutadas) ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.

A figura 1.1 mostra que o padrão geral dos dados não é linear na escala da função de ligação, assim utilizamos GAM. Possibilitamos até 1 intercepto por modelo neutro por sítio de amostragem (MN|SiteCode); para modelos com variáveis de dispersão continua não foi possível um smoother por Sítio de Amostragem e Modelo Neutro

1.1 GLMM binomial

1.1.1 Modelo Cheio e Auditoria

##                           dAICc   df  weight
## d/L_plot d/L_plot*MN|Site     0.0 22  1     
## k0 k0*MN|Site              7135.7 22  <0.001
## kf MN|Site                35074.5 123 <0.001
## d/L_plot MN|Site          60043.1 15  <0.001
## k0 MN|Site                61464.5 15  <0.001
## kf 1|Site                 65651.6 121 <0.001
## d/L_plot 1|Site           91136.5 13  <0.001
## k0 1|Site                 92878.5 13  <0.001

O único modelo plausível é aquele com função de ligação logito, a variável d/L_plot e a estrutura aleatória MN * d/L_plot|SiteCode.

Figura 1.2 Resíduos Quantílicos do modelo cheio para n_nRef ~ pdMN + p^2dMN + (d*MN|SiteCode)

NOTA: não entendo o motivo de não rodar os resíduos quantílicos quando compilo o Rmarkdown:

mensagem de erro:

Error in if (family$family %in% c(“binomial”,“poisson”,“quasibinomial”, : argumento tem comprimento zero Calls: … withCallingHandlers -> withVisible -> eval -> eval -> simulateResiduals Além disso: warining message: In checkModel(fittedModel) : DHARMa: fittedModel not in class of supported models. Absolutely no guarantee that this will work

1.1.2 Modelo cheio: predições

  cbind(n_nRef, 100 -
n_nRef)
Predictors Odds Ratios CI p
(Intercept) 37.07 13.35 – 102.96 <0.001
p.z 1.42 0.64 – 3.15 0.384
d_Lplot.z 0.59 0.34 – 1.03 0.063
MN [EI] 0.17 0.07 – 0.39 <0.001
p.z^2 0.30 0.13 – 0.69 0.005
p.z * d_Lplot.z 1.48 0.97 – 2.26 0.069
p.z * MN [EI] 0.34 0.17 – 0.66 0.001
d_Lplot.z * MN [EI] 0.97 0.43 – 2.19 0.944
d_Lplot.z * p.z^2 0.96 0.61 – 1.51 0.861
MN [EI] * p.z^2 1.04 0.52 – 2.09 0.915
(p.z * d_Lplot.z) * MN
[EI] 		0.50 0.27 – 0.94 0.033
(d_Lplot.z * MN [EI]) *
p.z^2 		1.07 0.55 – 2.07 0.848
Random Effects
σ2 3.29
τ00 SiteCode 11.53
τ11 SiteCode.d_Lplot.z 3.03
τ11 SiteCode.MNEI 7.96
τ11 SiteCode.d_Lplot.z:MNEI 6.94
ρ01 -0.02
-0.45
0.39
ICC 0.84
N SiteCode 103
Observations 4120
Marginal R2 / Conditional R2 0.105 / 0.861

Figura 1.3 Predito e observado por SiteCode

Figura 1.4 Probabilidade predito pelo modelo mais plausível

1.2 GLMM binomail subset:MN==EE

##              dAICc   df weight
## d + d^2|Site     0.0 15 1     
## d|Site        2588.4 12 <0.001
## 1|Site       20757.2 10 <0.001

O modelo mais plausível considera um termo linear e quadrático para d na estrutura aletória. Segue seleção de variáveis.

Figura 1.6 Resíduos Quantílicos do glmm cheio para MNEE

Seleção de Variáveis

# dados
df_md <- filter(df_resultados,MN=="EE")
# modelo global
global_md <- glmer(cbind(n_nRef,100-n_nRef) ~
                     ( p.z + I(p.z^2) ) * ( d.z + I(d.z^2) ) +
                     (d.z + I(d.z^2)|SiteCode),
                   family = "binomial",data=df_md,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),na.action = "na.fail")
dredge_glmmEE <- dredge(global_md,trace = FALSE,evaluate=TRUE)
subset(dredge_glmmEE,delta<7)
## Global model call: glmer(formula = cbind(n_nRef, 100 - n_nRef) ~ (p.z + I(p.z^2)) * 
##     (d.z + I(d.z^2)) + (d.z + I(d.z^2) | SiteCode), data = df_md, 
##     family = "binomial", control = glmerControl(optimizer = "bobyqa", 
##         optCtrl = list(maxfun = 1e+05)), na.action = "na.fail")
## ---
## Model selection table 
##     (Int)     d.z   d.z^2    p.z  p.z^2 d.z:p.z d.z:I(p.z^2) I(d.z^2):p.z
## 30  3.856 -0.3820         0.4117 -1.179  0.8548                          
## 32  3.884 -0.3320 0.07890 0.4111 -1.182  0.8581                          
## 62  3.861 -0.4017         0.4087 -1.183  0.8656     0.019820             
## 96  3.885 -0.3384 0.07500 0.3676 -1.184  0.7767                   -0.1387
## 64  3.888 -0.3493 0.07873 0.4084 -1.186  0.8674     0.017350             
## 160 3.889 -0.3312 0.09039 0.4094 -1.187  0.8608                          
## 224 3.955 -0.3349 0.21910 0.3254 -1.252  0.7652                   -0.2256
## 128 3.887 -0.3466 0.07493 0.3664 -1.186  0.7811     0.008296      -0.1385
## 22  2.678 -0.3783         1.1020         0.8500                          
## 192 3.891 -0.3459 0.08725 0.4076 -1.188  0.8671     0.014840             
## 256 3.943 -0.2402 0.24550 0.3323 -1.241  0.7112    -0.094730      -0.2430
## 24  2.706 -0.3298 0.07660 1.1030         0.8531                          
##     d.z^2):I(p.z^2 df    logLik    AICc delta weight
## 30                 11 -6409.908 12841.9  0.00  0.369
## 32                 12 -6409.775 12843.7  1.76  0.153
## 62                 12 -6409.906 12844.0  2.02  0.134
## 96                 13 -6409.357 12844.9  2.95  0.084
## 64                 13 -6409.773 12845.7  3.78  0.056
## 160      -0.011000 13 -6409.774 12845.7  3.78  0.056
## 224      -0.142600 14 -6409.047 12846.3  4.35  0.042
## 128                14 -6409.358 12846.9  4.98  0.031
## 22                 10 -6413.507 12847.1  5.18  0.028
## 192      -0.007982 14 -6409.777 12847.8  5.81  0.020
## 256      -0.168200 15 -6408.992 12848.2  6.27  0.016
## 24                 11 -6413.383 12848.9  6.95  0.011
## Models ranked by AICc(x) 
## Random terms (all models): 
## 'd.z + I(d.z^2) | SiteCode'

Tabela 1.2.2

Predito para cada Sítio

Para avaliar o ajuste do modelo ao conjunto de dados utilizo funções do pacote MuMIn (REF) que permite calcular a predição do modelo médio para o conjunto original de dados.

Figura 1.7 Predições por sítio de amostragem do modelo médio calculado a partir do modelo global ; pontos: número de SADs neutras refutadas para uma bateria de simulações com um determinada distância média de dispersão; eixo x = distância média de dispersão / largura da área de amostragem; a linha é a probabilidade de não refutar uma SAD neutra segundo a predição média do conjunto de sub-modelos dentro do intervalo de plausibilidade de 7 (Burnham et al 2011)

Predito para novo conjunto de dados

Para avaliar o predito e intervalo de confiança de 95% pelo modelo médio utilizo funções do pacote AICcmodavg (REF).

Figura 1.8 Efeito predito de d e p na Probabilidade de não refutar segundo o modelo médio para MN==EE. No painel superior a predição média, nos paineis inferiores o intervalo de confiança.

COISAS PARA FAZER FIG 18: i) melhorar a figura ii) comparar com o mesmo gráfico para o predito e IC bootsrap para o modelo cheio [SUGESTÃO POR EMAIL PI]

1.3 n_nRef subset:MN==EI

##              dAICc   df weight
## d + d^2|Site     0.0 15 1     
## d|Site        7920.4 12 <0.001
## 1|Site       63552.4 10 <0.001

O modelo mais plausível considera um termo linear e quadrático para d na estrutura aletória. Segue seleção de variáveis.

Figura 1.9 Resíduos Quantílicos do glmm cheio mais plausível para MNEI

# dados
df_md <- filter(df_resultados,MN=="EI")
# modelo global
global_md <- glmer(cbind(n_nRef,100-n_nRef) ~
                     ( p.z + I(p.z^2) ) * ( d.z + I(d.z^2) ) +
                     (d.z + I(d.z^2)|SiteCode),
                   family = "binomial",data=df_md,
                   control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=100000)),na.action = "na.fail")
dredge_glmmEI <- dredge(global_md,trace = FALSE,evaluate=TRUE)
subset(dredge_glmmEI,delta<7)
## Global model call: glmer(formula = cbind(n_nRef, 100 - n_nRef) ~ (p.z + I(p.z^2)) * 
##     (d.z + I(d.z^2)) + (d.z + I(d.z^2) | SiteCode), data = df_md, 
##     family = "binomial", control = glmerControl(optimizer = "bobyqa", 
##         optCtrl = list(maxfun = 1e+05)), na.action = "na.fail")
## ---
## Model selection table 
##     (Int)    d.z  d.z^2       p.z   p.z^2 d.z:p.z d.z:I(p.z^2)
## 79  2.735        -2.147  0.052260 -0.9115                     
## 80  2.935 0.5378 -1.967  0.056100 -0.9083                     
## 208 3.151 0.5308 -2.186 -0.073260 -1.1270                     
## 207 2.952        -2.344 -0.068360 -1.1250                     
## 96  2.923 0.5179 -1.980  0.194700 -0.9040  0.3705             
## 112 3.154 0.8376 -1.981 -0.009226 -1.1350              -0.3201
## 224 3.151 0.5320 -2.184  0.063220 -1.1270  0.3650             
## 240 3.214 0.6970 -2.134 -0.074240 -1.1910              -0.1693
## 128 3.070 0.7209 -1.980  0.109200 -1.0510  0.2517      -0.2028
## 256 3.109 0.4198 -2.219  0.088050 -1.0840  0.4312       0.1148
##     I(d.z^2):p.z d.z^2):I(p.z^2 df    logLik    AICc delta weight
## 79       -1.1980                11 -6288.847 12599.8  0.00  0.180
## 80       -1.1980                12 -6287.921 12600.0  0.17  0.165
## 208      -1.0710         0.2145 13 -6287.193 12600.6  0.74  0.124
## 207      -1.0830         0.2025 12 -6288.288 12600.7  0.91  0.114
## 96       -1.0740                13 -6287.311 12600.8  0.98  0.110
## 112      -1.1340                13 -6287.345 12600.9  1.04  0.107
## 224      -0.9524         0.2140 14 -6286.677 12601.6  1.74  0.075
## 240      -1.0710         0.1598 14 -6287.097 12602.4  2.58  0.050
## 128      -1.0740                14 -6287.169 12602.5  2.72  0.046
## 256      -0.9315         0.2507 15 -6286.649 12603.5  3.71  0.028
## Models ranked by AICc(x) 
## Random terms (all models): 
## 'd.z + I(d.z^2) | SiteCode'

Tabela 1.3.2

Predito para cada Sítio

Para avaliar o ajuste do modelo ao conjunto de dados utilizo funções do pacote MuMIn (REF) que permite calcular a predição do modelo médio para o conjunto original de dados.

Figura 1.10 Predito por SiteCode a partir do modelo médio para MNEI.

Figura 1.11 Efeito predito de d e p na Probabilidade de não refutar segundo o modelo médio para MN==EE. No painel superior a predição média, nos paineis inferiores o intervalo de confiança.

diff_S = (S_obs - S_MN)/S_obs

Padrão Geral

Figura 2.1 diff_S = (S_obs - S_MN)/S_obs + (-1) * min(diff_S) + 0.01, inclinação_MNEE ~ 0 para todo k.

Parece que existem alguns outliers. De maneira geral MNEE apresenta boa congruência com os dados enquanto o comportamento de MNEI depende da dispersão: para k->100% o desvio aumento para ps extremos (bimodal/quadratica), apresentando boa congruência para proporção intermediárias de habitat, mas sempre subestimando S_obs. Com o aumento da capacidade de dispersão as estimavas de S se tornam mais próximas ao observado para p>0.5, superestimando o observado para valores elevados de p e k>0.5 o que leva a reduzir a porção de p em que MNEI faz uma boa aproximação.

U - taxa de especiação necessária para obter a riqueza observada no equilíbrio

Padrões Gerais

Figura 3.1 Padrões gerais de U_med: ~ k (% de propágulos até a vizinhança imediata da planta progenitora); ~ Stotal (riqueza observada na área amostral). E o padrão empírico encontrado nos dados Stotal ~ p

GAMM

Como os modelos com estrutura aleatória mais complexa levam muito tempo para estimar, irei começar com a estrutura aleatória mais simples para avaliar qual a melhor variável de dispersão e se for necessário refazer a seleção da estrutura aleatória.

Tabela de Seleção do Modelo Mais plausível:

##                    dAICc  df               weight
## d log gamma           0.0 121.264568083601 1     
## 1-k log gamma       174.1 125.474107334272 <0.001
## k.f log gamma       227.4 154.787008201268 <0.001
## d log normal        394.5 124.471358521882 <0.001
## d inverse gamma     766.2 127.68541117751  <0.001
## 1-k log normal      855.7 127.94024990325  <0.001
## k.f log normal      889.4 170.667471702588 <0.001
## 1-k inverse gamma   983.2 124.456458934862 <0.001
## k.f inverse gamma  1002.8 140.80259373167  <0.001
## d inverse normal   1136.4 121.957626137592 <0.001
## k.f inverse normal 1345.5 165.620770162204 <0.001
## 1-k inverse normal 1354.2 122.442084920132 <0.001
## 1-k id normal      2163.0 125.186614559034 <0.001
## d id normal        2189.1 123.301338427342 <0.001
## k.f id normal      2242.8 160.04673943848  <0.001

Segue gráficos diagnosticos do modelo mais plausível

Figura 3.2 Gráficos Diagnóstico do modelo mais plausível

O modelo apresenta boa congruência com os dados, porém parece existir outliers no entando.

Figura 3.3 Métricas de influência dos dados (broom::augment())

Exclui os dois sítios com valores de grande influência do conjunto de dados:

Figura 3.4 Grafico diagnostico do modelo mais plausivel sem os outliers.

Figura 4.5 Efeitos Parciais do modelo mais plausível

Parâmetros dos Modelos Neutros

d / L_plot

padrões gerais

Figura 5.1.1 d/Lplot ~ k

Descrição Estatística

## [1] "exp(coef):"
## (Intercept)       k0.95        k0.9       k0.85        k0.8       k0.75 
##   0.6715478   1.4153327   1.7360273   2.0081002   2.2669913   2.5268653 
##        k0.7       k0.65        k0.6       k0.55        k0.5       k0.45 
##   2.7936350   3.0726860   3.3755022   3.7066527   4.0756037   4.4953275 
##        k0.4       k0.35        k0.3       k0.25        k0.2       k0.15 
##   4.9791588   5.5515868   6.2545191   7.1436593   8.3360091  10.0630874 
##        k0.1       k0.05 
##  12.9390064  19.3670608

Figura 5.1.2 Gráficos Diagnósticos do modelo gam(d ~ k + s(DA,by=k), family=Gamma(link=“log”)). Exponenciação dos coeficientes parametricos estimados pelo modelo.

m’

padrões gerais

Figura 5.2.1 m’ ~ p (~k). A variável m’ é uma função de p, d e I(1/J)

Descrição Estatística

##        dAICc  df               weight
## model2    0.0 158.703955395811 1     
## model3  488.0 43               <0.001
## model1 7043.7 113.502207953068 <0.001

Figura 5.2.2 Diagnostico do modelo mais plausivel

theta

padrões gerais

Figura 5.3.1 Padrão Geral de theta considerando escala log

Descrição Estatística

Apendice: código extra

n_nRef: GAMM

Comparação de Modelos Cheios

l_md <- vector("list",6)
names(l_md) <- c("k 1|SiteCode", "k MN|SiteCode",
                 "k_1 1|SiteCode","k_1 MN|SiteCode",
                 "d/L 1|SiteCode","d/L MN|SiteCode")
# muitos parametros para poucos dados "k_1 k_1 * MN|SiteCode","d_Lplot d_Lplot * MN|SiteCode")
# k
l_md[[1]] <- gam(cbind(n_nRef,100-n_nRef) ~ k * MN + s(p.z,k,by=MN,bs="fs") +
                   s(SiteCode,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
l_md[[2]] <- gam(cbind(n_nRef,100-n_nRef) ~ k * MN + s(p.z,k,by=MN,bs="fs") +
                   s(SiteCode,by=MN,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
# k_1
l_md[[3]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN + s(p.z,by=MN,bs="tp") + s(k_1.z,by=MN,bs="tp") + ti(p.z,k_1.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
l_md[[4]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN +  s(p.z,by=MN,bs="tp") + s(k_1.z,by=MN,bs="tp") + ti(p.z,k_1.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,by=MN,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
# d/L
l_md[[5]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN + s(p.z,by=MN,bs="tp") + s(d_Lplot.z,by=MN,bs="tp") + ti(p.z,d_Lplot.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
l_md[[6]] <- gam(cbind(n_nRef,100-n_nRef) ~ MN +  s(p.z,by=MN,bs="tp") + s(d_Lplot.z,by=MN,bs="tp") + ti(p.z,d_Lplot.z,by=MN,bs=c("tp","tp")) +
                   s(SiteCode,by=MN,bs="re"),
                   family = "binomial",
                   data = df_resultados, method = "REML")
AICctab(l_md,weights=TRUE)
##                 dAICc   df               weight
## k MN|SiteCode       0.0 567.902990229433 1     
## k_1 MN|SiteCode 10456.4 255.87941997766  <0.001
## d/L MN|SiteCode 14414.2 255.878481862603 <0.001
## k 1|SiteCode    29130.6 480.710534787585 <0.001
## k_1 1|SiteCode  38331.6 162.939249908097 <0.001
## d/L 1|SiteCode  43420.3 162.955609731983 <0.001
# save(l_md,file="~/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/Support Information/md_nRef.Rdata")
# load("~/Documentos/Doutorado/artigo_mestrado/1A.P_MOVER/Support Information/md_nRef.Rdata")

O único modelo plausível considera a variável k (fator) e estrutura aleatória MN|SiteCode.

##       dAICc df               weight
## by=MN   0.0 567.902990229433 1     
## by=k   47.0 497.198353791482 <0.001

Diagnostico do Modelo Cheio

Figura 1.2 appraise(gam( cbind(n_nRef,100-n_nRef) ~ k*MN + s(p.z,MN,by=k), binomial(logit) )

O modelo mais plausível não apresenta boa congruência com os dados: a distribuição binomial não parece se adequar bem aos dados e o predito e observado podem ser muito divergentes.

Figura 1.3 Probabilidade de não refutar uma SAD neutra para cada modelo estatístico (título dos paineis); na esquerda o modelo mais plausível, na esquerda o segundo modelo plausível (>1000 deltaAIC). Os valores médios de cada modelo neutro apresentam semelhante entre os dois gráficos, contudo o modelo mais plausível apresenta maior flexibilidade na descrição dos modelos.

Figura 1.4 Observado e predito pelo modelo mais plausível. Eixo x = prop. de cobertura vegetal na paisagem; Eixo y: número de SADs não refutadas. Em branco o padrão geral estimado e em vermelho o erro padrão.

o modelo não apresenta bom ajuste aos dados:
i) MNEE: de maneira geral a tendência global parece estar de acordo com o observado. Contudo, para k<0.5, a estimativa parece ser muito influenciada por pontos extremos.
ii) MNEI: para todo k, a tendência global não parece apresentar boa congruência com o observado principalmente quanto aos valores médios observados, e.g., para k>0.50 a tendência global parece subestimar o observado. Para k<0.5 há um ponto que parece ser de grande influência (k=0.25 p->0.5).

Figura 1.5 residuals.gam(l_md[[“k MN|SiteCode”]]) ~ p (~k*MN).

Os resíduos do tipo ‘deviance’ não indicam erro na estimativa quanto à tendência global; o quê para mim está errado principalmente para MNEI.

Melhoria do modelo mais plausível: CONSTRUÇÃO

Opções

Comparação metodo para estimar sp e classe de bs

smoother type

##    dAICc  df               weight
## cr    0.0 576.311752135652 1     
## cc 2939.1 543.980281194929 <0.001
## tp 3290.1 567.902990229433 <0.001
## ps 3506.6 559.143708020173 <0.001
## cp 6113.9 566.50730949499  <0.001
##         dAICc df               weight
## GACV.Cp   0.0 582.038923755289 0.5   
## GCV.Cp    0.0 582.038923755289 0.5   
## REML     38.0 576.311752135652 <0.001
## ML       38.4 576.287674101362 <0.001

Figura 1.6 Predito pelo gam com method=GACV.Cp e bs=cr. O intervalo de confiança foi removido por apresentar valores negativos.

Apesar das mudanças na configuração o modelo continua estimando MNEI errado.

Resumo dos achados da comparação:

  1. GCV.Cp e GACV.Cp são os métodos mais plausíveis para o smoothing type padrão (“tp”) e para o mais plausível (“cr”).
  2. Porém estes métodos de estimativa do parâmetro de smoothing apresentam ajustes estranhos aos padrões observados.
  3. O smoother type “mais plausível”cr" ainda apresenta o mesmo problema do anterior: MNEI não está fazendo bom ajuste.
  1. Como não consegui configurar os métodos GCV.Cp e GACV.Cp para um padrão mais adequado sigo com o terceiro mais plausível (“REML”, o padrão).
  2. Vou explorar ajustar modelos separados para cada conjunto de dados

funções de ligação para distribuição binomial

##         dAICc  df               weight
## cloglog    0.0 572.962297266108 1     
## logit   2249.8 576.311752135652 <0.001
## probit  3640.8 575.395171500865 <0.001

Figura 1.7 Gráficos Diagnósticos do modelo mais plausível binomial(link=cloglog)

A qualidade do modelo não melhorou.

Figura 1.8 Predição pelo modelo com função de ligação “cloglog”

Para MNEI o problema persiste. Além disso para MNEE parece que o ajuste piorou.

Avaliação do conjunto de dados separadamente

MNEE
##                        dAICc  df               weight
## k 1|SiteCode logit        0.0 261.404710989706 1     
## k 1|SiteCode probit     663.6 261.341234025727 <0.001
## k 1|SiteCode cloglog   1254.7 256.432441274999 <0.001
## d/L 1|SiteCode logit   7212.7 127.452558167134 <0.001
## d/L 1|SiteCode probit  7212.7 127.452558167134 <0.001
## d/L 1|SiteCode cloglog 7212.7 127.452558167134 <0.001
## k_1|SiteCode cloglog   7865.4 127.890778649655 <0.001
## k_1|SiteCode probit    8539.6 127.940128196497 <0.001
## k_1|SiteCode logit     8650.6 127.878762400437 <0.001

Figura 1.9 Gráficos Diagnósticos do modelo mais plausível para o conjunto de dados MN:EE

Figura 1.10 Predito para o conjunto de dados MNEE

Parece que alguns pontos tem grande grande influência nos dados. Sigo comparação de diferentes tipos de smoothers

##    dAICc  df               weight
## tp    0.0 270.017818977703 0.907 
## ts    4.6 255.449435050652 0.093 
## cr   84.5 261.404710989706 <0.001
## cs   94.1 254.247921741065 <0.001
## ps  340.7 259.020966845455 <0.001
## cc  947.0 240.615840486311 <0.001
## cp 1230.2 251.280447452311 <0.001

Figura 1.11 Graficos Diagnosticos para MNEE e bs=tp

Figura 1.12 Predito para MNEE e bs=“tp”

Ainda observo valores de grande influência.

MNEI
##                        dAICc   df               weight
## k 1|SiteCode cloglog       0.0 255.938425616834 1     
## k 1|SiteCode logit      3522.0 265.551522242843 <0.001
## k 1|SiteCode probit     4252.3 272.189573651799 <0.001
## k_1|SiteCode cloglog    6003.1 127.972940662994 <0.001
## k_1|SiteCode logit     10150.9 127.977463144225 <0.001
## k_1|SiteCode probit    10727.2 127.984432980891 <0.001
## d/L 1|SiteCode logit   15681.2 127.9410891041   <0.001
## d/L 1|SiteCode probit  15681.2 127.9410891041   <0.001
## d/L 1|SiteCode cloglog 15681.2 127.9410891041   <0.001

Figura 1.13 Gráficos Diagnósticos do modelo mais plausível para o conjunto de dados MN:EI

## modelo mais plausível
df_pred.k <- expand.grid(SiteCode=levels(df_resultados$SiteCode)[1],
                          k=levels(df_resultados$k),
                          p.z=seq(min(df_resultados$p.z),max(df_resultados$p.z),length=103)) %>% 
  mutate(p = p.z * sd(df_resultados$p) + mean(df_resultados$p))
df_pred.k <- cbind(df_pred.k,
                   predict(l_md.EI[["k 1|SiteCode logit"]],
                           df_pred.k,type="response",se.fit=TRUE))
df_pred.k[,c("fit","se.fit")] <- df_pred.k[,c("fit","se.fit")]*100

ggplot(filter(df_resultados,MN=="EI"),aes(x=p,y=n_nRef)) + 
  geom_point(alpha=0.5) + facet_wrap(~k,ncol=4) +
  # k
  geom_ribbon(aes(ymin=I(fit - 2*se.fit), ymax=I(fit + 2*se.fit), x=p),
              data=df_pred.k,
              alpha=0.3,
              inherit.aes=FALSE,
              color="red") +
  geom_line(aes(y=fit),color="white", data=df_pred.k) +
  ylab("Pr(not ref. SAD)") + xlab("Proportion of Tree Cover")

Figura 1.14 Predito para o conjunto de dados MNEI

O problema de estimativa errada persiste. Vou avaliar outros tipos de smoothers.

##    dAICc  df               weight
## cr    0.0 255.938425616834 1     
## cs   25.3 254.87197636468  <0.001
## cc 2237.8 256.227920344716 <0.001
## ps 2384.7 249.567127726108 <0.001
## tp 2698.6 250.18188320226  <0.001
## ts 2701.8 243.207853913116 <0.001
## cp 4728.6 261.426662661009 <0.001

O smoother type “cr”, modelo utilizado até o momento, é o mais plausível. Vou avaliar se a remoção de alguns pontos de grande influência podem melhorar o ajuste.

–>

–>

–> –>

–> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –>

–>

–>

–>

–> –>

–> –> –> –> –>

–>

–>

–>

–> –> –>

–>

–> –> –> –> –> –> –> –> –>

–>

–>

–> –> –> –> –> –> –> –>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –>

–> –> –> –>

–>

–> –> –> –>

–>

–> –> –>

–>

–> –> –> –> –> –> –> –> –>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –>

–> –> –> –>

–>

–> –> –>

–>

–> –>

–>

–>

–>

–> –> –>

–>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –> –>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –>

–>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –>

–>

–> –> –>

–> –> –> –> –> –>

–>

–> –> –> –> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–>

–>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –>

–>

–>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –> –>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –>

–>

–>

–> –>

–> –>

–> –>

–> –> –> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–> –> –> –> –>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–> –> –> –> –> –> –>

–>

–> –> –>

–>

–>

–>

–>

–> –> –>

–>

–>

–> –> –> –> –> –> –> –> –> –> –> –>

–>

–>

–>

–> –> –> –> –> –>

–>

–>

–> –> –> –> –> –> –> –>

–> –> –> –>

–>

–>

–> –> –> –>

–>

–>

–> –> –> –>

–>

–> –> –> –> –>

–>

–> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –>

–> –>

–> –>

–> –> –> –> –> –> –> –>

–> –>

–> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –>

–> –>

–> –> –> –> –> –> –> –>

–> –>

–> –>

–> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –> –>

–> –>

–> –>